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ABSTRACT 



Aims. We investigate the stability, nonlinear development and equilibrium structure of vortices in a background shearing 
Keplerian flow 

Methods. We make use of high-resolution global two-dimensional compressible hydrodynamic simulations. We introduce 
the concept of nonlinear adjustment to describe the transition of unbalanced vortical fields to a long-lived configuration. 
Results. We discuss the conditions under which vortical perturbations evolve into long-lived persistent structures and 
we describe the properties of these equilibrium vortices. The properties of equilibrium vortices appear to be independent 
from the initial conditions and depend only on the local disk parameters. In particular we find that the ratio of the 
vortex size to the local disk scale height increases with the decrease of the sound speed, reaching values well above the 
unity. The process of spiral density wave generation by the vortex, discussed in our previous work, appear to maintain 
its efficiency also at nonlinear amplitudes and we observe the formation of spiral shocks attached to the vortex. The 
shocks may have important consequences on the long term vortex evolution and possibly on the global disk dynamics. 
Conclusions. Our study strengthens the arguments in favor of anticyclonic vortices as the candidates for the promotion 
of planetary formation. Hydrodynamic shocks that are an intrinsic property of persistent vortices in compressible 
Keplerian flows are an important contributor to the overall balance. These shocks support vortices against viscous 
dissipation by generating local potential vorticity and should be responsible for the eventual fate of the persistent 
anticyclonic vortices. Numerical codes have be able to resolve shock waves to describe the vortex dynamics correctly. 

Key words, accretion, accretion disks - planet formation - vortices - hydrodynamics - methods: numerical 



1. Introduction 

Recent advances in the understanding of vortex behav- 
ior in differentially rotating flows are mainly associated 
with the study of p roto planetary disk dyn a mics, since 
Ivon Weiszackerl (| 19441 ) and I Adams fc Watkinsl (|l995l ), who 
suggested that vortices can promote the formation of plan- 
ctcsimals. Indeed, it has been shown that vortices, if sus- 
tained long enough, lead to particle aggregation in their 
core ( see e.g. [C havanis 2000; dc la Fuentc Marco s fe Bared 



2001; Johansen et al. 2004; Klahr & Bodenheimer 200 



and to the formation of protoplanets. 

The vortex scenario for planetary formation encounters 
an apparent obstacle: any structure in a Keplerian disk is 
subject to a strong shearing that may eventually lead to its 
decay. The only mechanism for sustaining a stable vortex 
in such flows is nonlinearity. Hence, vortices that may start 
the process of planetary formation should exceed a critical 
threshold in their amplitude. Direct numerical simulations 
are therefore an important tool in these studies. 

Previous studies on vortex s tability, performed main! 
by numerical simulations (see 



200C; Davis 



Godonfe Livid Il999b 



ability, performed mainly 
Barre fc Sommcria l!995t 



Klahr fc Bodenheimer 



iDavis et al 



2003; Klahr 



2004] iBarranco fc Marcudl2005HUmurhan fc Regevll2004h 
have shown that only anticyclonic vortices can be non- 
linearly stable in Keplerian flows and can promote dust 
trapping and the creation of protoplanets. Similar indica- 



tions come from hydrodynamical simulations of accretion 
disks with initially imposed random vort icity fields, where 
again only anticyclonic y ortices survive (|Shen et all [2006: 
I Johnson fc Gammid [20051 ). 

Although different numerical simulations have already 
shown that anticyclonic vortices can survive in compressible 
Keplerian disks, most of the stable structures are observed 
on scales smaller than the scale height of the disk. This calls 
for three-dimensional investigations be fore solid conclusions 
can b e drawn. On the other hand, IBarranco fc Marcus! 
(2005) have found indications of a 3D vortex instability 
and have shown that only off mid-plane vortices with a very 
specific vertical structure are stable in 3D. The preferred 
location for dust capture and planetesimals are, however, 
the high density regions in the midplane of the disk, where 
vortices are unstable. A solution to this problem could be 
the existence of vortices with a horizontal extent larger than 
the disk thickness. Those should in fact behave as 2D struc- 
tures therefore being stable and promote mass accumula- 
tion. However, presently, there are dou bts on the possibility 
of the existence of such vortice s (see iJohnson fc Gammid 
120051: IB arranco fc Marcus! [20051) . They are thought to de- 
cay to a smaller size, of the order of the disk scale height, 
due to the supersonic velocity circulation. 

In the present paper, we study the nonlinear dynamics 
of relatively large scale (exceeding the disk half thickness) 
vortices in keplerian flows by global compressible Simula- 
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tions. In this situation, compressibility is an important fac- 
tor that cannot be neglected because the Keplerian pro- 
file, setting a strong radial velocity shear, co uples neces- 
sarily vortex and wave mode perturbation (see iBodo et al.l 
|2005| ). Our aim is to study in detail equilibrium struc- 
ture of nonlinear long lived vortices. In our study we as- 
sume the initial vortical perturbation as given and we 
investigate how it evolves into an equilibrium configura- 
tion and the properties of this configuration. It is not our 
aim here to study how such large amplitude perturbation 
was f o rmed. Recent studies (see e.g. iKlahr fc Bodenheimer I 
l2003t iPetersen et al.l l2007f ) have pointed out how baro- 
clinic effects can lead to a growth of vortical perturbations, 
but that a radial temperature variation is required. 
However, here, for studying the vortex structure, 
we assume, for simplicity, a disk with initially uni- 
form temperature. The adequacy of this approxi- 
mation must be tested through future studies We 
will focus our analysis on the possible existence of vortices 
with horizontal extent larger than the disk scale height. An 
important process for the vortex dynamics in compressible 
flows is t he linear generatio n of spiral density waves, de- 
scribed in lBodo et al.l (|2005l ). In the present paper we fur- 
ther investigate this process by verifying its efficiency dur- 
ing the nonlinear stages and by analyzing its relevance for 
the vortex structure. In this respect, an interesting point is 
that equilibrium vortex structures work as a mill that pro- 
cesses the shear flow energy into coherent wave emission. 
The nonlinear development of spiral-density waves leads to 
the formation of spiral shocks with a steady pattern. These 
shocks may increase the stability of anticyclonic vortices by 
slowing down their decay and may have also global effect 
on the disk. 

In $2] we describe the numerical code, setup and initial 
conditions used in the simulations. In we present the 
results of our numerical simulations and describe the non- 
linear adjustment to an equilibrium configuration. We con- 
sider the effect of the vortex initial amplitude and size on its 
further evolution and analyse the timescales of the adjust- 
ment process. In 2] we discuss the stability and structure of 
the developed vortex structures. In JjHwe study the process 
of wave emission induced by the vortex and the subsequent 
shock development. A study on the effects of numerical res- 
olution is given in fj5J The paper is finally summarized in 

m 



where v r and v,p — L ( / ) /(pr) are the radial and azimuthal 
velocities, p is the thermal pressure and g r = — 1/r 2 is the 
gravitational acceleration. The total energy E of the fluid 
is expressed as the sum of the kinetic and internal energy 
via the ideal equation of state: 



E = 



T - 1 



(5) 



where T = 5/3 is assumed. Dissipative effects are only 
those related to the use of a mesh of finite width Equations 
([T])-((4]) are sol ved using the hydrody namics module of the 
PLUTO code (|Mignone et al.l l2007f ). PLUTO evolves the 
conservative equations in a Godunov-type fashion by first 
interpolating volume averages inside each cell and then 
solving a Riemann problem at each interface to compute 
the fluxes. For the present work, we choose second order 
slope-limited interpolati on and the two-stage R unge Kutta 
time stepping scheme of iGottlieb fc Shul (|1998f ). An accu- 
rate Riemann solver based o n the two-shock approximation 
()Colella &: Woodwardlll984[ ) is used to compute numerical 
fluxes. 

Integration of the equations is made cons iderably faste r 
by taking advantage of the FARGO scheme (Massct 2000) , 
available in the PLUTO code. This scheme allows larger 
timesteps than the standard integration by removing the 
average azimuthal velocity from the Courant condition, 
which severely limits the time step because of the fast or- 
bital motion at the inner boundary. Besides the increased 
computational efficiency, this technique considerably re- 
duces the amount of numerical diffusion in the solution. 

Our computational domain covers the region < 4> < 
2ir, R m < r < i? ou t, where R ln and i? ut are the inner- 
most and outermost boundaries in the radial direction. The 
grid is constructed by first assigning the number of 
uniform zones in the </> direction, whence A(f> = 2tt/N^,. 
The mesh spacing in the radial direction is then chosen 
to yield approximately equal radial and azimuthal lengths 
Ar^ = riAip on each ring i = 1, • ■ • , N r of the disk. The 
number of radial zones N r is determined by the condition 
ri + Ari/2 = r^+i — Ar^+i/2 (grid continuity) at each node, 
giving a first (non-integer) estimate 



n; = log 



/log 



(6) 



2. Numerical code and setup 

We solve the equations of inviscid compressible gas dynam- 
ics in two dimensional polar coordinates (r, 0) . From the 
observer's frame, let p, m r , and E be, respectively, the 
fluid density, radial, angular momentum and total energy 
density; their conservation is expressed by 



dp 1 d{rpv r ) ^ 1 d{pv 4> ) 
dt 



d{pv r ) 1 d{rpv 2 T ) 



dr r d(f> 

1 d(pv r v^) dp 

dt r dr r d(f> dr 
dL^ 1 d(rL<pv r ) 1 d^^v^ + rp) 
dt r dr r d<p 

dE_ + l d[r{E+p)v r ] + l d[{E + p)v^\ 
dt r dr r d<p 



0. 


(1) 


p4 4 

r 


" P9r , (2) 


o, 


(3) 


m r g r . 


(4) 



Finally, we obtain N r by rounding N* up to the nearest 
integer which requires i? ou t to be properly adjusted by in- 
verting Eq. ([6]) with N* replaced by N r . Hence, the radial 
resolution is set up by the azimuthal spacing and by the 
innermost and outermost boundaries. 

The initial condition, at t = 0, consists of a Keplerian 
disk with uniform density and pressure plus a velocity per- 
turbation described in the following subsection. 

Boundary conditions are set as follows: at the outer- 
most radius free outflow is permitted, whereas at the in- 
nermost radius variables are kept constant to their initial 
value. Periodic boundary conditions hold in the azimuthal 
domain. 

For our simulations we used four main numerical setups, 
with low, medium, high and very high resolutions (called 
respectively L, M, H and VH), described in detail in Table 

HI 
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Table 1. Resolution, radial domain and number of revo- 
lutions at r = 1 for the different numerical setups. The 
domain length in the angular direction is € [0, 2tt]. 





Resolution 


Domain 


Time 




(JV* x N r ) 


[-Rin, -Rout] 


(Revol.) 


L 


1500 x 550 


[0.2, 2] 


80 


M 


2000 x 622 


[0.4, 4] 


200 


H 


4000 x 1466 


[0.2, 2] 


80 


VH 


8000 x 1559 


[0.5, 1.7] 


80 



2.1. Initial vortex 

Initial conditions for numerical simulations are composed 
by the sum of the Keplerian flow with = r -1 / 2 and 
a vortex perturbation, that can have different geometries, 
size and amplitude. The initial vortex configuration most 
widely used in our simulations is the following: 



<(0) = ± 



■ exp 



~ x Q ) 2 (y - y ) 2 ) 



q 2 a 2 



w y(0) = Ttq{x - x ) exp 
p'(0) = 0. 



(x - x ) 2 (y - y ) 2 ) 



q 2 a 2 



(7) 



(8) 



(9) 



Here e defines the amplitude of the initial perturbation and 
its sign determines the vortex polarity (positive in the case 
of anticyclonic vortex). The parameters a and q describe, 
respectively, the size (in the radial direction) and aspect 
ratio of an elliptic vortex. A circular vortex corresponds 
to q = 1, and q > 1 refers to a vortex elongated in the 
azimuthal direction. We choose the location of the vortex 
(xq, 2/0 ) so that it corresponds to the radial location Vq = 1. 
With these choices, our unit of length is the radius tq at 
which the vortex is located and our unit of time is the 
inverse of the Keplerian angular velocity f2Ke P (?"o) a * the 
same location. From a physical point of view, the vortex 
amplitude e is related to the maximum vorticity f2jvf (at 
the vortex center), Q,m = 2e 

A second type of initial vortex with algebraic distribu- 
tion of potential vorticity is also used, 



<(0) = ± 



e(y - yo) 



l 



(x-x ) 2 , (y-y ) 2 ) 



? ? 

q*a z 



v'(0) = Teq(x-x o ) 1 + 



(x-x ) 2 , {y-yo) 2 ) 



q 2 a 2 



(10) 



.(11) 



This initial vortex is used for the sake of comparison to see 
how different initial potential vorticity fields with similar 
amplitude and scale but different geometry behave in the 
nonlinear regime. 

3. Vortex evolution and nonlinear adjustment 

3.1. The simulations 

The evolution of the perturbations depends on three main 
parameters: the amplitude and size of the perturbation (re- 
spectively e and a) and the sound speed in the disk, c s . Our 
first aim is to determine the region, in the parameter space 
described by e, a and c s , in which the evolution of the initial 



Table 2. Table of parameters used for the computations. 
From the left to right, the columns give the sound speed 
(c s ), the vortex initial size (a), amplitude (e), ellipticity 
parameter q, geometry (namely the shape of the initial vor- 
ticity distribution: exp and alg stand for the initial configu- 
rations given by Eqns. l7l8l and llOU"Tl respectively) together 
with its polarity Pol (— for anticyclonic vortices and + for 
cyclonic ones). The numerical setup (NS, given in Table 
[lj and an identification label are given in the rightmost 
columns. 
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c. 
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Vjreum/ jrui 




T nhpl 
ijdUcl 


U.x 


U. 1 


U. 1 


i 
i 


exp/ - 


T 


A 1 

ill 






n o 
u.z 


i 

i 


exp/ - 


T 

1^ 








n q 
U . o 


i 

i 


exp/ - 


T 

Lj 








0.5 


x 


exp/ - 


H 


A4 






n s 

U . o 


i 

i 


exp/ - 


T 






0.3 


0.5 


1 


PYTl / - 


L 


A6 




0.4 


0.5 


1 


OXT) / - 


L 


A7 


yJ.XJO 


u.o 
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exp/ - 


M 


Bl 


n ni 

U.Ul 


n ni 

U.Ul 


u . 


i 
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exp/ - 


V 1 1 








i 

i . 


i 

i 


exp/ - 


u 
n 






0.02 


0.5 


1 


OXT) / - 


H 


C3 




0.05 


0.5 


1,2 


exp,alg/ - 


VH 


C4 




0.1 


0.1 


1 


exp/ - 


H 


C5 






0.2 


1 


exp/ - 


H 


C6 






0.22 


1 


exp/ - 


H 


C7 






0.25 


1 


exp/ - 


H 


C8 






0.28 


1 


exp/ - 


H 


C9 






0.3 


1 


exp/ - 


H 


C10 






0.5 


1 


exp/ - 


H 


Cll 




0.2 


0.2 


1 


exp/ - 


M 


C12 






0.5 


1 


exp/ -,+ 


H 


C13 


0.005 


0.02 


0.5 


1 


exp/ - 


H 


Dl 




0.025 


0.5 


1 


exp/ - 


H 


D2 


0.002 


0.02 


0.5 


1 


exp/ - 


H 


El 


0.001 


0.01 


0.5 


1 


exp/ - 


H 


E2 




0.02 


0.5 


1 


exp/ - 


H 


E3 




0.05 


0.5 


1 


exp/ - 


H 


E4 




0.1 


0.5 


1 


exp/ - 


H 


E5 



perturbation leads to a stable, long-lived equilibrium config- 
uration. As we shall see, before reaching this final state the 
system undergoes, in the course of several disk revolutions, 
a transition phase that we call nonlinear adjustment. Our 
second aim is to provide a detailed description of the equi- 
librium vortex configuration. With these purposes we have 
performed runs with different values of the three param- 
eters and different numerical setups, using Eqs. [7][8] with 
q = 1 (circular vortices). In particular, by increasing the 
value of the sound speed (c s = 0.001,0.01,0.1), we have 
explored the behavior of the vortex changing a at fixed e 
and changing e at fixed a. Additional calculations have been 
performed for the purpose of exploring in more detail par- 
ticular regions of the parameter space. For instance, addi- 
tional values of c s have been used for a better understanding 
of the scaling behaviors of some of the vortex properties. 
Moreover, in order to investigate whether the general be- 
havior is changed by varying the shape and structure of the 
initial perturbation, we have performed computations with 
different values of the ellipticity parameter q and of the (al- 
gebraic) potential vorticity distribution (see Sect, above). A 
complete list of all the simulations and related parameters 
are described in Table [2] 
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3.2. Nonlinear amplitude thresholds 

It is well known that vortices with small amplitudes are 
sheared away by the linear deformation induced by the 
background flow. Nonlinear effect can counteract the shear- 
ing deformation only for vortices with higher amplitude. In 
this respect, we found two threshold values for e defining 
the character of the evolution of the initially imposed vor- 
tex. The first threshold value between linear and nonlinear 
behaviors is: 

e*=0.1. (12) 

For e < e* the vortex behaves linearly (Bod o et alj 12005). 
Its dynamics is characterized by the geometrical stretch- 
ing of the initial configuration, due to the strong radial 
Keplerian velocity shear. The linear vortex is then stretched 
to fill the entire (2ir) azimuthal domain, within the first 
1.5 local disk revolutions. It then experiences further decay 
down to the dissipation scales, where it is damped. 

When the vortex amplitude exceeds the first nonlinear 
threshold, a two-stage process occurs. First, the vortex is 
sheared into a narrow vortex layer, which then undergoes 
local instabilities. We then observe the formation of small- 
scale weak anticyclonic vortices at different azimuthal loca- 
tions. Hence, the initial vortex is transformed into a chain 
of small scale vortices that tend to spread and occupy the 
entire azimuthal domain. The typical character of the evo- 
lution of potential vorticity for this case is shown in Fig. 

ED 

Increasing the initial vortex amplitude, we found a sec- 
ond nonlinear threshold: 

e** = 0.25. (13) 

Vortices with e > e** experience direct adjustment from 
the initial to the final persistent structure, i.e. a strong an- 
ticyclonic vortex. This process may be followed in Fig. [2l 
where an initially unbalanced vortex, within 4 revolutions, 
undergoes direct adjustment to the final nonlinear config- 
uration. The energy excess of the initial state is radiated 
away in the form of spiral-density waves and shocks. 

The final equilibrium configuration appears to be a non- 
linear attractor reached by the system if the initial ampli- 
tude exceeds e** and if the initial spatial scale falls in a 
range discussed in the next subsection. Indeed, the same 
nonlinear state is developed from all initial vortices satis- 
fying these conditions, independently from the details of 
the initial potential vorticity distribution (exponential or 
algebraic, circular or elliptic). 

3.3. Initial vortex size 

The result of nonlinear adjustment strongly depends also on 
the initial vortex size. In order to quantify this dependence, 
we have carried out computations with vortices of different 
initial size a. The size of the vortex can be compared with 
the length-scale 

H = c s /n Kcp , (14) 

which, in presence of the vertical component of gravity, de- 
scribes the scale height of a thin Keplerian disk. From a 
physical point of view, a 2D analysis can give meaningful 
results only for vortices with sizes larger than H. This will 
be our case of interest. For c s initially constant all over 



the disk (as it is in our case), H is a function of the radial 
position and, at the vortex location, one has H — c s (since 

n Kep (i) = i). . .. 

Our numerical results show that vortices undergoing 
direct nonlinear adjustment have initial size in a range 
ao < a < a max , where the limits ao and a max depend 
on the local disk parameters. During the nonlinear adjust- 
ment the vortex decreases its size, radiating the excess en- 
ergy through spiral density waves and shocks, and reaches 
the equilibrium configuration with size ao that can be still 
larger than the disk scale height. The final equilibrium con- 
figuration reached by these vortices appears to be inde- 
pendent from the details of the initial state and seems to 
represent a nonlinear attractor, whose characteristics will 
be described in detail in the next Section. In cases C1-C13 
(i.e., c s = 1CP 2 ), the value range of a evolving into the same 
nonlinear configuration is 2H < a < 10H. For comparison, 
in cases E2-E5 (c s = 10~ 3 ) this range is 12H < a < 2QH. 

When the spatial-scale of the initial vortex exceeds 
QWx, the evolution is quite complex. We observe a radial 
transfer of potential vorticity in both directions (see Fig 
[3]) , caused by the combined action of shocks (radiated from 
the initially imposed supersonic vortex) and of flow curva- 
ture (inducing Rossby wave variations). The initial radial 
transfer of potential vorticity is accompanied by a shearing 
process leading to a decrease of potential vorticity local- 
ization and to a subsequent fragmentation. The resulting 
picture tends to become even more complicated when sec- 
ondary vortices, formed at different radii and with different 
sizes, start to interact. As a result, we see multiple vortices 
with typical sizes equal to or smaller than H leaning to 
further decay in size and amplitude. 

When the size of the vortex is smaller than ao and 
the initial vortex amplitude exceeds the second nonlinear 
threshold, we observe nonlinear adjustment to a final con- 
figuration with sizes similar to the initial values. In this 
case, developed vortices typically have sizes smaller than 
ao. However, if the initial amplitude largely exceeds the 
second nonlinear threshold, we have found that the vor- 
tex can increase its size and reach the nonlinear attractor 
discussed above. We have not studied this case (a < ao) 
in detail since the typical size of vortices developed during 
the adjustment falls around H or smaller (unless the initial 
amplitude is quite high), indicating the importance of 3D 
consideration. 

Fig. Q] shows a sketch diagram illustrating the character 
of vortex evolution depending on its initial amplitude (e) 
and size (a). 

In the following, we will focus our attention on the cases 
that reach the equilibrium configuration. For them we dis- 
cuss, in the following subsection, the time-scales of nonlin- 
ear adjustment and, in the following 21 the final equilib- 
rium structure. 

3.4. Time-scales of the adjustment 

The characteristic time-scale needed for radiating the en- 
ergy excess and reaching the equilibrium configuration 
strongly depends on the vortex size. As we shall see in the 
next Section, the equilibrium vortex configuration has a 
characteristic size that does not depend on the initial size 
of the vortex. This means that bigger initial vortices need 
to emit more energy to adjust to the final small-scale con- 
figuration. Indeed, we have seen that a different number 
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P.Vo'ticity; t=1 rev. Density; t=1 rev. 




-0.50 -0.40 -0.30 -0.20 -0.10 0.00 0.10 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 



Fig. 1. Dynamics of an initially imposed small-scale anticyclonic vortex (exponential distribution of vorticity, q = 1, 
a = 0.1) with amplitude exceeding the first nonlinear threshold parameter: e = 0.2. The system parameters correspond 
to case C6 (see Table 2). The first column shows the evolution of the perturbed potential vorticity and the second column 
the evolution of density at times 1, 5 and 10 (in units of the local disk revolution period). The third column displays 
zooms at the vortex location at the same times. The vortical perturbation is initially stretched into a narrow vortex sheet 
which then undergoes local nonlinear instabilities producing a chain of small-scales anticyclonic vortices that tend to fill 
the entire azimuthal domain. See electronic version of the paper for color figures. 

of revolutions were needed to develop a stable long-lived structure starting from vortices of different size: 

a = 0.05 — * 3 revolutions, 
a = 0.1 — » 4 revolutions, 
a = 0.2 — > 6 revolutions 1 . 
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P.Vorticity; t=1 rev. 



Density; t 



y; t= I rev. 



-1 
P.Vorticity; t=8 rev. 



Density; t=8 rev. 



P.Vorticity; t=20 rev. 



-1 
Density; t=20 rev. 



■ ■ 






-1 o 



.00 -0.78 -0.55 -0.55 -G.1C 0.12 0.55 -0.60 -0.28 0.03 0.35 0.67 0.98 1.30 



Fig. 2. Dynamics of an initially imposed small-scale anticyclonic vortex (exponential distribution of vorticity, q = 1, 
a = 0.1) with amplitude exceeding the second nonlinear threshold parameter: e = 0.5. The system parameters correspond 
to case Cll (see Table 2). The first column shows the evolution of the perturbed potential vorticity and the second column 
the evolution of density at times 1,8 and 20 revolutions (in units of the local disk revolution period). The initially imposed 
vortex with high amplitude emits the excess energy in form of density-spiral waves and within 4 revolutions adjusts to 
the single nonlinear configuration which survives further without significant variation throughout the entire period of 
simulation. See electronic version of the paper for color figures. 



(runs C4,C11 and C13). The geometry of the potential vor- terns observed in prev ious numerical simulations (see e.g. 
ticity adjustment can be seen in Fig. [51 It matches the pat- iGodon fc Livioll 1999a ). However, our high- resolution simu- 
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Fig. 3. Potential vorticity evolution in the case of an anticyclonic vortex with exponential distribution of vorticity, q = 1, 
a = 0.2 and amplitude exceeding the second nonlinear threshold parameter: e = 0.5. The system parameters correspond 
to case C13 (see Table 2). The background Keplerian vorticity is subtracted. The dark spots reveal the anticyclonic 
patches. The radial transport of potential vorticity occurs on the initial stages, when the vortex is oversized and is 
subject to curvature radiation. At later times the background shearing stretches the vorticity and developed vortices 
only move azimuthally, while occasionally interacting and decaying. 



lations reveal new features of this process. The adjustments 
of potential vorticity and density proceed in different ways 
and on different timescales (see Fig. [6j). When the potential 
vorticity has almost reached its equilibrium configuration, 
the density distribution still undergoes significant struc- 
tural variations. Initially it develops a double core vortex 
configuration that slowly tends to change into a one core 
stable vortex. For an initial value of a = 0.1, the poten- 
tial vorticity adjusts in 4 revolution, while density needs 15 
revolutions for reaching the final configuration (case Cll). 
Moreover, the density adjustment (that is, the time needed 
for the merging of the two density cores) shows a strong 
dependence on the intrinsic numerical viscosity. Indeed, in- 
creasing the resolution, this settlement time increases. The 
differences between vorticity and density adjustments may 
be a consequence of the existence of spiral shocks gener- 
ated by the interaction of the nonlinear vortex with the 
background shear flow, as we will discuss in Sj5l 

4. Vortex stability and structure 

One of the main goals of the present study is to describe the 
stability and structure of long-lived vortices in Keplerian 
disks. For this purpose, we selected cases undergoing direct 
adjustment to a single vortex, and we followed their long 
time behavior. Fig [7] shows radial profiles of potential vor- 
ticity and density at the center of a vortex for case Cll, 
integrated for a total of 200 revolutions. After 4-6 rotations 
the potential vorticity configuration shows a well defined 
center. The density maximum matches the vorticity center 
after approximately 15 revolutions. We observe a profound 
persistent structure throughout 200 revolutions. Of partic- 
ular interest is the behavior of the mass accumulated by the 
anticyclonic velocity circulation, shown in Fig. [71 where we 
observe a slow but steady increase of density at the vortex 
center. The same figure shows the steepening of the vortic- 
ity gradient during the adjustment process. 

By analyzing vortices with different initial configura- 
tions, we have found that the shape of the developed, per- 
sistent vortex depends only on the disk sound speed c s , and 
not on the characteristics of the initial vortex (provided its 
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Fig. 4. Sketch diagram illustrating the character of the evo- 
lution of an initially imposed anticyclonic vortex. In ab- 
scissa we have the amplitude (e) and in ordinate the size (a) 
of the initial anticyclonic vortex. Two nonlinear thresholds 
in amplitude (0.01 and 0.25) control the fate of the initial 
vortex, as it can decay linearly, evolve into a vortex chain 
or undergo direct adjustment to a single persistent configu- 
ration. The initial vortex size determines if the vortex will 
decay due to radial vorticity transport (oversized vortices), 
viscosity, or undergo direct adjustment. Horizontal dotted 
line shows the size of the equilibrium persistent configura- 
tion which can be interpreted as the nonlinear attractor. 



size a and amplitude e fall in the ranges discussed in the 
previous Section). We can characterize the final state by its 
size ao (in the radial direction) and ellipticity parameter q. 
The vortex size is measured defining as vortex the region 
that has potential vorticity larger than 10% of the max- 
imum value at the vortex center. The following equation 
give s a scaling law for op as fu nction of the sound speed c s 
(see iBarranco fc Marcusj[2005l) : 

a = f(c s )H, (15) 

where the function f(c s ) (i.e. the ratio ao/H) is shown 
in Fig. O Note, that although our vortex scales with the 
sound speed (and the disk height scale), we have found 
vortices with maximal size larger than the ones described 
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Fig. 5. Scaling law for the equilibrium vortex size: f(c s ) 
(= clq/H) as a function of sound speed c s . The radial extent 
of the vortex is used. 

in lBarranco fe Marcus! (|2005f ). The plot shows that the ra- 
tio ao/H increases for decreasing the sound speed, so the 
condition for a 2D vortex behavior (vortex size larger than 
the disk scale height), is verified to a better extent at lower 
sound speeds. The final ellipticity q, on the other hand, ap- 
pears to be almost independent from c s and has a typical 
value ~ 5. 

We confirm that the initial vortex configuration, for giv- 
in g direct adjustmen t, has to satisfy the condition derived 
in lDavis et alj (|2000f ): the velocity gradient in the center of 
the initial vortical field must exceed the Keplerian gradient. 
Moreover, we observe that the velocity gradient of the per- 
turbation tends to decrease during its evolution and reaches 
the Keplerian value in the equilibrium configuration. 

Fig. [5] shows the decay curves of the maximum of po- 
tential vorticity at the center of the vortex for different 
resolutions. As expected, for finer meshes, the vortex de- 
cay is slower because of the reduced numerical dissipa- 
tion. However, apart from this consideration, interesting 
effects can be noticed. After a rapid variation within the 
adjustment time (first 5 revolutions in present case), one 
can observe the formation of a stable persistent vortex. 
Interestingly, the vortex seems to be temporarily able to 
oppose viscous dissipation, exhibiting for some time an in- 
crease of the maximum potential vorticity. This effect im- 
plies production of potential vorticity. On the other hand, 
potential vorticity is a nonlinearly conserved quantity in 
barotropic flows. Thus, the source of the coherent gener- 
ation of potential vorticity necessary to support or even 
enhance the vortex are the spiral shock waves that will be 
discussed in They seem to be able to overcome the (nu- 
merical) dissipation at early times, but not at later times. 
The situation is quite complex, one of the possible reasons 
for this behavior could be a change in the internal vortex 
structure, like the transition from the double to the single 
density core. 

4.1. Cyclonic vortices 

The nonlinear dynamics of cyclonic vortices shows signif- 
icant differences from anticyclonic ones. In the limit of 
small amplitudes, cyclonic vortices decay on the shearing 
timescale, related to the geometrical stretching induced by 
the background differential rotation. For higher amplitudes 
(exceeding the first nonlinear threshold e*) they exhibit a 



rapid decay on time-scales shorter than the shearing time- 
scales. It then appears that nonlinear forces accelerate the 
linear decay of cyclonic vortices. 

5. Waves and shocks 

Vortices in shear flows generate spira l density waves by 
a line ar mechanism first described by Chagelis hvili et all 
(|l997l ). and further invest igated numerically in Keplerian 
flows bv lBodo et alj (|2005[ ). In the present context we study 
the characteristics of this process when nonlinear forces 
are also in action. As clearly seen from Figs [TJ [5] and [51 
the dynamics of potential vorticity is accompanied by wave 
emission, not only during the transition time, but also af- 
terwards, when the nonlinear self-sustained vortex is fully 
developed. This enforce the fact that nonlinearity does not 
suppress the wave emission. Moreover, after the adjust- 
ment, waves emitted by the coherent vortex structure ap- 
pear to develop into spiral shocks. 

A regular structure of spiral shock waves has been found 
in protop lanetary disk simulations with an embedded pro - 
toplanet ijTanigawa fc Watanabdl2002t iKoller et al.ll2003l ). 
Traces of shock generation by vortices may be found in 
Uohnson fe Gammid (|2005l) . In our simulations we verify 
the existence of a steady pattern of spiral shocks produced 
by a single vortex. Moreover, shock waves appear to be an 
inherent property of vortices in sheared compressible flows, 
but they can be observed only at fairly high resolutions us- 
ing shock capturing schemes (limits on resolution will be 
discussed in fQ. Our high resolution calculations allow to 
study these shock waves in great detail. In Fig [8] we show 
the distributions of potential vorticity, density, temperature 
and local Mach number for the vortex and the attached 
shock waves. One can distinctly recognize a wave-crest of 
the density-spiral wave developing into a double shock con- 
figuration, with the shock ahead (behind) of the vortex fac- 
ing the outward (inward) region. A couple of much weaker 
shocks, parallel to the strong ones, appear to be present 
although they remain barely visible. These shocks strongly 
affect the density structure of the developed vortex configu- 
ration, resulting in a splitting of the vortex core. Eventually, 
however, the shearing background leads to the merging of 
the cores, but the shocks persist. 

Spiral shocks induced by a wake of planets are believed, 
in s ome situations, to be respons ible for planet migration 
(see iPapaloizou fe Terqueml 120061 and references therein) . 
In our computations no radial variation of the vortex posi- 
tion has been observed. 

As seen in previous studies, spiral shocks affect dust 
accretion rates on the vortex core and thus promote the 
formation of planetesimal. In this sense, they increase the 
importance of anticyclonic vortices in planetary formation. 

The presence of shocks has consequences for the vortex 
evolution. Nevertheless, the final fate of these structures 
cannot be easily foreseen and requires much longer simu- 
lations. We can here only sketch some possible scenarios. 
One possibility is the exhaustion of matter in the vortex 
bearing ring and th e formation of isolated planetesimal. 
ISchafer et alj (|2004l ) have argued that spiral shocks may 
lead to the gap formation. On the other hand, shocks heat 
the ring at the radius where the vortex is sustained (see 
Fig. [5]) , which in turn may trigger the linear Rossby wave 
instabili ty due to the unusua l entropy gradient in the disk 
matter (|Lovelace et al.l Il999f) . The instability will induce 
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Fig. 6. Evolution of potential vorticity during the direct nonlinear adjustment. The initial conditions correspond to a 
small-scale anticyclonic vortex with e = 0.5, a = 0.05 and q = 1. The system parameters correspond to case C4 (see 
Table 2). Potential vorticity of the initially imposed vortical perturbations undergoes adjustment to the final equilibrium 
state within the first 3 revolutions. 
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Fig. 7. Evolution of potential vorticity and density for the anticyclonic vortex with e = 0.5, a = 0.1, q — 1 (setup Cll). 
Radial cuts are taken across the geometrical center of the vortex (corresponding to a minimum in potential vorticity, for 
anticyclonic vortices). The figures show a gradual increase of density in the vortex core. Viscous forces damp oscillations 
around the vortex. At higher resolution this effect occurs on longer time-scales. After the initial adjustment, the vortex 
persists in the flow for long times and its dynamics is only a subject to viscous damping. 



radial mixing and possibly the destruction of the coher- 
ent vortices. On the other hand, iLubow fc Qgilvid {l998) 
have shown that spiral shocks can be themselves unstable 
in three dimensions. Hence, as we said, longer and possibly 
three-dimensional simulations can clarify this issue. 

6. Resolution study and numerical requirements 

Our direct numerical experience shows that a number of 
stringent numerical requirements may be indispensable for 
observing the correct dynamics of small-scale vortices in 
two-dimensional compressible Keplerian flows. One neces- 
sary requirement is the need to resolve small size perturba- 
tions. 

Vortices resolved on 8 computational zones are persis- 
tent and, although affected by enhanced numerical viscos- 
ity and damping, still survived more than 20 revolutions. 
Otherwise, at least 24 grid points are necessary to cor- 
rectly follow the nonlinear adjustment. In Fig [TO] we show 
an example of such a situation, when numerical setup M 
(2000x733, 6 grid points over vortex) failed to describe the 
direct nonlinear adjustment of the initially imposed small- 
scale vortex, but displayed development of a dipolar or even 
tripolar vortex configurations. Direct nonlinear adjustment 
to the single long-lived structure is resolved only under con- 
siderably higher resolution (setup VH, 8000x1559, 24 grid 
points over the vortex) thus setting overall restriction on 
the resolution of our code. 



Notwithstanding the fact that the purpose of the sim- 
ulations is the study of small-scale vortices, global simula- 
tions are needed: vortices interact globally and change the 
density, temperature and potential vorticity in the whole 
ring at the radius where they are situated. Hence, to get 
a correct dynamical picture, we need true global simula- 
tions with a domain covering 2-7T in azimuthal direction. 
Calculations on smaller domains in <j> (tt/2) have lead to 
drastically different results, since the two vortices at the 
same radius start to interact before the nonlinear adjust- 
ment occurs. 

Another requirement on the numerical code is that it 
should be able to capture shocks, an important ingredient 
of the overall balance of the sustained structures in differen- 
tially rotating flows. In our case the PLUTO code employs 
solvers that accurately resolve the shock dynamics. This 
requirement does not play in favor of spectral codes for 
compressible simulations of the Keplerian disk dynamics. 



7. Summary 

Direct numerical simulations of the nonlinear dynamics of 
coherent vortices in 2D compressible protoplanetary disks 
with Keplerian differential rotation have been presented. 

We have shown the possible existence of anticyclonic 
vortices with sizes exceeding the Keplerian disk height 
scales. We have followed the evolution of such vortices for 
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Fig. 8. Potential vorticity, density, temperature and local Mach number of the anticyclonic vortex after 10 revolutions 
in the disk. Here the initial vortex has e = 0.5, a — 0.05, q — 1 (case C4). Narrow spiral rays beaming out of the 
central vortex in the potential vorticity panel reveal spiral shocks attached to the outer radius of the vortex. Density 
perturbations show that compression induced by these shocks are responsible for the double core appearance of the 
vortex. The temperature plot shows the heating induced by shock dissipation. Heating is stronger in the vicinity of the 
vortex leading to a consistent temperature increase of the global azimuthal ring that bears the vortex. See electronic 
version of the paper for color figures. 



200 local revolutions, showing their persistence and stabil- 
ity. 

We have found that the development of long-lived non- 
linear anticyclonic vortex configurations occurs only when 
the amplitudes of the initially imposed vortex perturbations 
exceed some threshold value. We have distinguished two 
nonlinear threshold values, that define the fate of the ini- 
tial vortex. At small amplitudes, the vortex behaves linearly 
and it decays due to the shearing deformation. When the 
amplitude exceeds the first threshold, the vortex stretching 
is followed by nonlinear instabilities, leading to the forma- 
tion of multiple small-scale anticyclonic vortices. At ampli- 
tudes higher then the second threshold, direct transition 
from the initially unbalanced to the final equilibrium con- 
figuration occurs. We have interpreted the latter process as 
the nonlinear vortex adjustment and studied the param- 
eters that can describe this process. We found that the 
vortex adjustment proceeds on two different time-scales. 
Shorter times are required for the settling the potential 
vorticity (3-6 revolutions) and longer time-scales are re- 
quired for the adjustment of density. Density adjustment 
time-scales depend on the amount of viscosity present in 
the flow. Thus, we observe a transient double-core vortex, 
that evolves into a single core structure. The density core 
splitting of the vortex is due to the spiral shock induced 
compression. Cores tend to merge at later times and the 
process is promoted by the intrinsic numerical viscosity. At 
lower viscosity, double core vortices survive for longer times. 



The density at the center of the vortex slowly increases and 
the potential vorticity gradient steepens. 

The structure of the developed long-lived vortex does 
not depends on the initial vortex configuration, provided 
it exceeds the second threshold amplitude and its size does 
not exceed a limiting value. In this sense we found a nonlin- 
ear attractor that is the final configuration of a wide range 
of initial vortical perturbations. We found a scaling law for 
the equilibrium vortex size. In particular we find that the 
ratio of the vortex size to the local disk scale height in- 
creases with the decrease of the sound speed. Therefore, 
at low sound speed the radial extent of equilibrium anticy- 
clonic vortices can significantly exceed the disk thickness. 
In this case vortex dynamics is intrinsically 2D and should 
be well modelled by the present simulations. Note that the 
sound speed required is somewhat lower than the typical 
values expected for protoplanetary disks (c s ~ 0.05 — 0.1). 
In the expected range the radial extent would be of the or- 
der of the scale height, while the azimuthal extent will be 
five times larger, in this case, to have a conclusive answer, 
3D calculations are needed. 

Vortices generate density-spiral waves that rapidly de- 
velop into shocks. As a result, a long-lived nonlincarly 
balanced vortex is accompanied by two spiral compress- 
ible shock waves facing both radial directions. The shock 
strength is maximal in the vicinity of the vortex outer edge, 
where the shock is attached to the main vortex circulation 
and heats all the annular disk region containing the vortex. 
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Higher resolution runs indicate that these shock waves are 
able to generate local potential vorticity and may support 
the vortex against viscous dissipation. 

We analyzed cyclonic vortices at nonlinear amplitudes. 
It seems that the linear decay due to the shearing deforma- 
tion is accelerated by nonlinear effects. 

We performed a resolution study and found the neces- 
sity of a global high-resolution approach to numerical sim- 
ulations of these processes. The fate of the vortex depends 
significantly on the dissipation properties of the code. So, 
in this sense the use of physical viscosity, as opposed to 
the intrinsic numerical viscosity present in our simulations, 
would improve understanding of the vortex behavior and 
damping. 

Our study contributes to the scenario of planetary for- 
mation inside the core of the long-lived vortices. We found 
that protoplanetary disks with lower sound speed can sus- 
tain vortices with higher ratio of vortex size to disk thick- 
ness and create a more favorable conditions for dust trap- 
ping and mass accumulation. In this context, we have also 
found a steady increase of density inside the nonlinearly 
balanced vortex, partly, due to the existence of persistent, 
steady, spiral shock waves that we showed to be an intrinsic 
property of stationary vortices with sizes exceeding the half 
thickness of the disk. 
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time (revolutions) 

Fig. 9. Plots of the maximum value of potential vorticity 
(st the vortex center) vs time, at three different resolutions. 
After the initial adjustment (vertical dotted line) vortex 
opposes viscous dissipation for some period and falls after- 
wards into the almost linear decay regime. At 2000 points 
in the azimuthal direction (12 grid points per vortex scale) 
vortex decays in about the 200 revolutions, while at 4000 
point resolution the estimate of the vortex disappearance is 
at 350 revolutions. At 8000 points the lifespan of the vortex 
is longer, but can not be confidently estimated using our 
data. 
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Fig. 10. Dynamics of an initially imposed small-scale an- 
ticy clonic vortex with q = 1, a = 0.01 and e = 0.5 at two 
different resolutions: (2000 x 733, 6 grid points over the vor- 
tex) on the left and (8000 x 1559, 24 grid points) on the 
right. The fate of the vortex adjustment strongly depends 
on the resolution. At lower resolution we see the develop- 
ment of a bipolar vortex. Fully resolved vortex development 
needed resolution as high as 8000 points covering global az- 
imuthal domain. 
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